Using scRNA-seq to identify markers of rice root cell differentiation
Here, I’ll work through the pre-processing and then analaysis of a publicly available dataset of 10x Genomics scRNA-seq data from rice roots. These data were prepared by the Wang Lab at the Institute of Plant Physiology and Ecology and downloaded from the NCBI bioproject PRJNA706435 and SRA run SRR13853440.
setwd("~/Desktop/scRNAseq-rice/")
library(DropletUtils)
library(tidyverse)
library(ggpointdensity)
library(scico)
library(scales)
library(irlba)
library(Matrix)
library(Seurat)
library(tibble)
library(patchwork)
# These funcs are copied from a kallisto/bustools tutorial made by
# the Patcher Lab (https://www.kallistobus.tools/tutorials/)
get_knee_df <- function(mat) {
total <- rank <- NULL
tibble(total = Matrix::colSums(mat), # total reads per droplet
rank = row_number(dplyr::desc(total))) %>% # ranks droplets in desc order by total
distinct() %>%
dplyr::filter(total > 0) %>%
arrange(rank)
}
get_inflection <- function(df, lower = 100) {
log_total <- log_rank <- total <- NULL
df_fit <- df %>%
dplyr::filter(total > lower) %>%
transmute(log_total = log10(total),
log_rank = log10(rank))
d1n <- diff(df_fit$log_total)/diff(df_fit$log_rank)
right.edge <- which.min(d1n)
10^(df_fit$log_total[right.edge])
}
knee_plot <- function(df, inflection) {
total <- rank_cutoff <- NULL
annot <- tibble(inflection = inflection,
rank_cutoff = max(df$rank[df$total > inflection]))
ggplot(df, aes(total, rank)) +
geom_path() +
geom_vline(aes(xintercept = inflection), data = annot, linetype = 2,
color = "gray40") +
geom_hline(aes(yintercept = rank_cutoff), data = annot, linetype = 2,
color = "gray40") +
geom_text(aes(inflection, rank_cutoff,
label = paste(rank_cutoff, "'cells'")),
data = annot, vjust = 1) +
scale_x_log10() +
scale_y_log10() +
labs(y = "Rank", x = "Total UMIs") +
annotation_logticks()
}We’ll start by loading the cell x gene count matrix that’s output from Kallisto/ Bustools. These files are in the repository along with the steps used to make them.
# load output from kallisto/bustools
# matrix where rows are genes, cols are cells
res_mat <- ReadMtx(mtx="cells_x_genes.mtx",
cells="cells_x_genes.barcodes.txt",
features="cells_x_genes.genes.txt",
feature.column = 1,
mtx.transpose = TRUE)
# Read in genes and gene names
tr2g <- read_tsv("t2g.txt", col_names = c("transcript", "gene", "gene_name",
"drop", "chromosome", "start", "stop", "strand"))
tr2g <- distinct(tr2g[, c("gene", "gene_name", "chromosome")])# Inspect library saturation
# Empty/near-empty droplets will have ~0 reads and ~0 observed genes
lib_sat <- tibble(nCount = colSums(res_mat), # total reads per droplet
nGene = colSums(res_mat > 0)) # total number of genes observed in a droplet
ggplot(lib_sat, aes(nCount, nGene)) +
geom_bin2d(bins = 50) +
scale_fill_distiller(palette = "Blues", direction = 1) +
scale_x_log10() + scale_y_log10() + annotation_logticks()
# It seems like the data from the SRA accession may be pre-filtered for quality
# Normally you would expect to see high density around 0,0Another way to assess library saturation is with a knee plot. On the y-axis is droplet/cell rank, determined by the total number of reads in a droplet (x-axis). The inflection point of the curve indicates a reasonable threshold for filtering droplets/cells by total read count.
knee_df <- get_knee_df(res_mat)
inflection <- get_inflection(knee_df)
# again, based on the shape of the curve this data seems to be pre-filtered
knee_plot(knee_df, inflection)
# Filter out low readcount cells (likely empty or lower-quality droplets)
res_mat <- res_mat[, colSums(res_mat) > 500]
res_mat <- res_mat[Matrix::rowSums(res_mat) > 0,]We also want to filter out droplets with a high percentage of reads mapped to either the mitochondria or chloroplast. Droplets where these reads are more abundant may indicate droplets in which the nuclear transcripts of lysed cells are more likely to have been lost/degraded, but organellar reads persisted because they were protected by an additional membrane layer.
# get a list of genes from organellar genomes
mt <- filter(tr2g, chromosome=="Mt")$gene %>%
intersect(., rownames(res_mat))
pt <- filter(tr2g, chromosome=="Pt")$gene %>%
intersect(., rownames(res_mat))
# create Seurat obj and calculate percent reads mitochondrial/chloroplast per cell
seu <- CreateSeuratObject(res_mat)
seu[["percent.mt"]] <- PercentageFeatureSet(seu, features=mt)
seu[["percent.pt"]] <- PercentageFeatureSet(seu, features=pt)
# this plot looks OK
# would look worse if low readcount cells hadn't already been filtered above
p.mt <- ggplot(seu@meta.data, aes(nCount_RNA, percent.mt)) +
geom_pointdensity() +
scale_color_scico(palette = "batlowW", direction = 1, end = 0.9) +
labs(x = "Total read count", y = "Percentage mitochondrial") +
theme_bw()
p.pt <- ggplot(seu@meta.data, aes(nCount_RNA, percent.pt)) +
geom_pointdensity() +
scale_color_scico(palette = "batlowW", direction = 1, end = 0.9) +
labs(x = "Total read count", y = "Percentage chloroplast") +
theme_bw()
p.mt | p.pt
# Remove cells with > 10% mitochondrial reads or 1.25% chloroplast
seu <- subset(seu, subset = percent.mt < 10 & percent.pt < 1.25)With the matrix now filtered, we normalize the read counts to reduce skew in read count distribution and stabilize variance across genes. There are plenty of ways to normalize, but here we use the Seurat default which is log normalization. Once normalized, we can look for genes with high variation in expression from cell to cell.
# Normalize data
seu <- NormalizeData(seu)
# identify a subset of genes that exhibit high cell-to-cell variation
seu <- FindVariableFeatures(seu)In RunPCA(), we use these variable features (rather than all >30K genes) as the basis for dimension reduction. We can also visualize the dimension loadings to see which genes explain the most variance on a given axis and make an elbow plot to show which axes explain the most variance across the dataset.
# run PCA
seu <- ScaleData(seu)
seu <- RunPCA(seu)
# visualize the dim loadings
p1 <- VizDimLoadings(seu, dims = 1:2, reduction = "pca") +
theme(axis.text = element_text(size=6))
p2 <- ElbowPlot(seu) + theme(axis.title.y = element_text(size=10))
(p1 | p2) + plot_layout(design="AAB") & theme(axis.text = element_text(size=9))
UMAP is a useful way of visualizing the dimension reduction done via PCA. Here, we use the first 15 axes of the PCA results to run UMAP (another dimension-reducing method) to visually represent those 15 axes in just 3 (or the number of n.components).
# run other dimension-reduction algorithms
seu <- RunUMAP(seu, dims = 1:15, n.neighbors = 100,
n.components = 3, n.epochs = 500)
DimPlot(seu, reduction = "umap", dims = c(1,2)) +
theme(legend.position = "none") +
labs(x = "UMAP 1", y = "UMAP 2")
We can then use these data to identify neighbors and clusters, i.e. cells that are closer together in the reduced-dimension space and therefore have similar overall gene expression patterns——and therefore indicate different cell types, treatments, etc.
seu <- FindNeighbors(seu, reduction="umap", dims=1:3,
n.trees = 500, k.param=100)
seu <- FindClusters(seu, resolution = 0.7, n.start=100, n.iter=1000)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 15711
## Number of edges: 1825821
##
## Running Louvain algorithm...
## Maximum modularity in 100 random starts: 0.9104
## Number of communities: 20
## Elapsed time: 41 secondsDimPlot(seu, reduction = "umap", dims=c(1,2)) +
labs(x = "UMAP 1", y = "UMAP 2")